library(sf) # vector manipulation
library(raster) # raster manipulation
library(fasterize) # "faster" raster
library(whitebox) # terrain analysis
library(dplyr)
library(rgdal)
library(tidyverse)
library(AOI)
library(osmdata) # OSM API
library(elevatr) # Elevation Web Tiles
library(gifski) # Gifski
FLooding Jaron
stage = height of the water in a channel streamflow = the rate of water flowing in a channel basin = an area in which all cells contribute to a common outlet (an area or ridge of land that separates waters flowing to different rivers, basins, or seas) flowpath = the linear path over which water flows (river) HAND = Height Above Nearest Drainage, “how high is a cell above its nearest river cell?”
#1-1.2 Collecting Data, Basin Boundary, Elevation Data
basin = read_sf("https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-11119750/basin/")
write_sf(basin, dsn = "/Users/xingxin/Github/geog176a-summer-2020-lab1/USGS-11119750.gpkg")
elev = elevatr::get_elev_raster(basin, z = 13) %>%
crop(basin) %>%
mask(basin)
elev_basin = elev * 3.281
writeRaster(elev_basin, filename = "/Users/xingxin/Github/geog176a-summer-2020-lab1/mission-creek-basin-elev.tif", overwrite = TRUE)
elev_raster = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/mission-creek-basin-elev.tif")
#1.3 Buildings and river-network data
bb = st_bbox(basin) %>%
st_as_sfc() %>%
st_transform(4326)
#Building
buildings = osmdata::opq(bb) %>%
add_osm_feature(key = "building") %>%
osmdata_sf()
#Railways
railway = opq(bb) %>%
add_osm_feature(key = 'railway', value = 'station' ) %>%
osmdata_sf()
#Streams
stream = osmdata::opq(bb) %>%
add_osm_feature(key = 'waterway', value = "stream") %>%
osmdata_sf()
#1.3 Buildings poly and river-network and railways
#buildings data
buildings_points = buildings$osm_lines %>%
st_intersection(basin) %>%
st_transform(crs(basin))
buildings_poly = buildings$osm_polygons %>%
st_intersection(basin) %>%
st_transform(crs(basin)) %>%
st_centroid()
#railways data
railways = railway$osm_points %>%
st_intersection(basin)
#streams data
streams = stream$osm_lines %>%
st_intersection(basin)
#2.1 Terrain Analysis
##Hillshade
##River-raster
wbt_hillshade("/Users/xingxin/Github/geog176a-summer-2020-lab1/mission-creek-basin-elev.tif","/Users/xingxin/Github/geog176a-summer-2020-lab1/wbt-hillside.tif")
## [1] "hillshade - Elapsed Time (excluding I/O): 0.15s"
hillshade = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/wbt-hillside.tif")
plot(hillshade, axes = FALSE, box = FALSE, col = gray.colors(256, alpha = 0.5), main = "Hillshade", legend = FALSE)
plot(streams$geometry, add = TRUE, col = "navy")
plot(basin$geometry, add = TRUE)
#2.2.1 Height Above Nearest Drainage - River Raster
river_raster = streams %>%
st_transform(5070) %>%
st_buffer(10) %>%
st_transform(crs(elev_raster))
stream_raster = fasterize::fasterize(river_raster, elev_basin)
writeRaster(stream_raster, filename = "/Users/xingxin/Github/geog176a-summer-2020-lab1/river-raster.tif", overwrite = TRUE)
##2.2.2 Height Above Nearest Drainage - Hydrologically Corrected Surface
streams_raster = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/river-raster.tif")
wbt_breach_depressions("/Users/xingxin/Github/geog176a-summer-2020-lab1/mission-creek-basin-elev.tif","/Users/xingxin/Github/geog176a-summer-2020-lab1/hydro-corr-surf-elev.tif")
## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.100s"
wbt_elevation_above_stream("/Users/xingxin/Github/geog176a-summer-2020-lab1/hydro-corr-surf-elev.tif","/Users/xingxin/Github/geog176a-summer-2020-lab1/river-raster.tif", "/Users/xingxin/Github/geog176a-summer-2020-lab1/wbt_elevation_above_stream.tif")
## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.28s"
##2.2.3 Height Above Nearest Drainage - HAND Raster
HAND = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/wbt_elevation_above_stream.tif")
HAND_raster = HAND + 3.69
streams_raster = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/river-raster.tif")
HAND_raster[streams_raster == 1] = 0
writeRaster(HAND_raster, filename = "/Users/xingxin/Github/geog176a-summer-2020-lab1/hand-raster.tif", overwrite = TRUE)
#3.1 2017 Impact Assessment
Floods = raster("/Users/xingxin/Github/geog176a-summer-2020-lab1/hand-raster.tif")
Flood_map = Floods
Flood_map[Flood_map > 10.02] = NA
#3.2 Plots
plot(hillshade, axes = FALSE, box = FALSE, col= gray.colors(256, alpha = .5), legend=FALSE)
plot(Flood_map, add=TRUE, col= rev(blues9), legend=FALSE)
plot(railways$geometry, add=TRUE, col= "green", cex=1, pch=16)
cols = ifelse(!is.na(raster::extract(Flood_map, buildings_poly)), "red", "black")
Flood_map[Flood_map >= 10.02] = NA
plot(hillshade, axes = FALSE, box = FALSE, col = gray.colors(256, alpha = 0.5), legend = FALSE, main = paste(sum(cols =="red"), "impacted building"," 10.02 foot stage"), cex = 0.5)
plot(Flood_map, add = TRUE, col = rev(blues9))
plot(buildings_poly, add = TRUE, col = cols, cex = .08, pch = 16)
plot(railways, add = TRUE, col = "green", cex = 1, pch = 16)
plot(basin$geometry, add = TRUE, border = "black")